Field effect on surface states in a doped Mott-Insulator thin film 
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Surface effects of a doped thin film made of a strongly correlated material are investigated both 
in the absence and presence of a perpendicular electric field. We use an inhomogeneous Gutzwiller 
approximation for a single band Hubbard model in order to describe correlation effects. For low 
doping, the bulk value of the quasiparticle weight is recovered exponentially deep into the slab, 
but with increasing doping, additional Friedel oscillations appear near the surface. We show that 
the inverse correlation length has a power-law dependence on the doping level. In the presence of 
an electrical field, considerable changes in the quasiparticle weight can be realized throughout the 
system. We observe a large difference (as large as five orders of magnitude) in the quasiparticle 
weight near the opposite sides of the slab. This effect can be significant in switching devices that 
use the surface states for transport. 

PACS numbers: 71.30.+h, 71.27. +a, 73.61.-r 
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I. INTRODUCTION 

The metal insulator transition (MIT) based on carrier 
doping of a Mott insulator has been investigated exper- 
imentally and theoretically^^. Recently, the formation 
of a superconducting phase was observed at the inter- 
face of a Mott and band insulator and the possible tun- 
ing of these transitions by an external electric field was 
reported 3 -. Moreover a three terminal setup was imple- 
mented by Son et al. who induced hole doping in a thin 
Mott insulator film from a doped band insulator through 
the application of a voltage difference between the drain 
and the gate terminals^. 

For the above class of phenomena inhomogeneities and 
proximity effects play an essential role. In order to deal 
with such systems one needs a theoretical model that 
is able to include correlation effects in heterostructures 
while not being too computationally expensive such that 
one has the possibility to consider large enough system 
sizes. This is crucial especially for the investigation of 
systems where surfaces and finite size effects are signifi- 
cant such as thin films made of strongly correlated ma- 
terials. The interface between a band insulator and a 
strongly correlated system has been studied theoretically 
with a two site dynamical mean field theory(DMFT)^ 
and the slave boson mean field theory(SBMFT}&. Such 
studies predict the formation of a two dimensional elec- 
tron gas at the interface which arises due to charge re- 
construction. Surface correlation effects were studied 
theoretically in half filled heterostructures modeled by 
a single band Hubbard model!. Also the penetration of 
metallic behavior into a Mott insulator was studied both 
within the Gutzwiller approximation and DMFT for the 
half filled case&2. Surface correlation effects of a doped 
semi-infinite Hubbard model were investigated within an 
embedded DMFT for both single band and multi-band 
system a 10 ' 11 . Within this method, due to numerical lim- 
itations, only few surface layers (up to 6) can be used in 
order to address site dependent correlation effects. When 
the correlation length is large, this method is not reliable 
any more. 



In order to describe position dependent electronic corre- 
lation effects in a slab geometry we employ the Gutzwiller 
approximation(GA). While GA works only for the metal- 
lic phase, it gives reliable information about the quasi- 
particle (QP) weight of electrons at different spatial lo- 
cations. For heterostructures, GA was found to be in 
good qualitative agreement with the more refined DMFT 
method for the half filled caseS. While GA and SBMFT 
are equivalent for zero temperature^, in two site DMFT, 
like GA, the bulk QP weight is governed by a simple 
power law and there is only a correction to U c when 
compared with the linearized DMF T 13 ' 14 . Generally, GA 
over-estimates the QP weight at all dopings but it is con- 
sidered to be accurate enough to describe low energy ex- 
citations and is routinely used for interpolations in com- 
bination with DMFT methods^. 

The aim of this paper is to investigate the spatial depen- 
dence of the charge density and the QP weight of a doped 
correlated slab and to understand the correlation effects 
in the presence of an external electric field. We predict 
significant changes in the QP weight throughout the sys- 
tem. This study is motivated by potential applications 
in nanoscale switching devices with spatial controllable 
conductivity through the application of an external elec- 
tric field. 

The outline of the paper is as follows: after a brief 
derivation of the saddle point equations for a slab geome- 
try(section II) the results for a doped correlated slab are 
presented in section III(A). Next the effect of an electric 
field is discussed in section III(B) and finally we present 
our conclusions in section IV. 



II. MODEL AND METHOD 

The simplest Hamiltonian that is able to capture the 
essential physics of strongly correlated systems is the sin- 
gle band Hubbard models, 
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where t%j are the hopping amplitudes, (ij) is summa- 
tion over nearest neigthbour sites and U is the Hubbard 
energy describing the Coulomb interaction between two 
electrons with opposite spin located on the same site. 
In the presence of an external electric field the model 
becomes^: 

H = Hu + ^Vih ia , (2) 

where Uj is the position dependent potential. In spite of 
the simple form of the Hubbard model, exact solutions 
exist only for d = 1 and d = oo^ i 17 i 18 and therefore we 
are forced to work with approximations. If one is only 
concerned about ground state properties or low energy 
excitations^, one of the choices is the Gutzwiller approx- 
imation(GA) which is the infinite dimension limit of the 
Gutzwiller wave function (GWF^i 2 ^. GWF is a many 
body wave function with an additional degree of freedom 
used to reduce the weight of higher energy configurations. 
In the single band Hubbard model these configurations 
are on-site double occupancies obtained when two elec- 
trons with opposite spin reside on the same site. The 
GWF is written as: 

\GWF)=l[P i \<t> ), (3) 

i 

where i is the lattice site index and the projector opera- 
tors are defined as Pj = g e) ie.i + # CT ,iS<r,i + 9»,iS&,i + 9d,id%- 
The operators e = (1 - ri la )(l - h is ), s a = h ia {\ - n is ) 
and d = fii^nis are local projectors of zero, singly and 
doubly occupied states, | v^o) is & noninteracting Fermi sea 
and consist of both spin up and spin down states and the 
g coefficients are variational parameters. The following 
local constraints have to be satisfied in order to remove 
the local contributions in the diagramatic expansion of 
various expectation value a 19 i 21 , 

(PjPi)o = 1, (4) 
(P}Pi4* c i*)o = kW)o- (5) 

Where (...)o represents the expectation value with re- 
spect to |</9o)- The explicit form of the above constraints 
is the following: 

g«r 2 (ei)o + ^.9icr 2 (s. 1C r)o + 9d/{di)o = 1, (6) 

a 

9ia 2 (sia)o + 9d/(di)o = (*W) - ( 7 ) 

In the limit of infinite dimensions the effect of the pro- 
jectors Pi requires the renormalization of the hopping 
amplitudes between different site s 18 i 21 . These renormal- 
ization factors can be written as: 




By substituting Eqs. (J5]) and ([7]) into Eq. ((5]) one arrives 
at an expression for y/q ia that is only a function of gd,i, 
(n ia ) and (h lS ) as: 



V(l - (ht)o + dj)((h ia ) - dj) + y/di((n is ) - dj) 
(nio-)o(l - (nicr)o) 



where di = g 2 j(n-jcr) 0(^3)0 is the probability of dou- 
ble occupancy that is calculated within \GWF) and 
(ni)o = (nia) + {hi S ) n . Moreover, in addition to Eq .([5]) 
the relation (rii CT ) gutzwiller = ("■i(r)o holds in the limit of 
infinite dimensions. By considering the above relations 
the total energy functional of an inhomogeneous system 
has the following form, 



(H)gwf = ^2 -^jv / gi7y / ^7(c| a -Cjcr)o + y^Vj (n icr ) Q, 

(ij),a i,<? 

+ yUg 2 dyl (h l(J ) (n lS ) 07 (9) 



The conditions (GWF\GWF) = {<p \<Pa) and {<p \<p ) = 
1 are used in the above relation, the first relation itself 
is a consequence of the infinite dimensional limit and the 
second relation is just the normalization condition for 

\<M- 

Away from half filling the problem of minimizing the 
energy functional is combersome because the renormal- 
ization factors, qi ai are functions of (fii a ) . Therefore 
it is impossible to simply vary the above energy func- 
tional with respect to (0o I ■ A possible approach, similar 
to DFT, is to start with an arbitrary value for (nj a .) 
and then to expand the energy functional as function of 
(ni a } up to linear order around the starting (fij cr ) . This 
allows us to vary the energy functional with respect to 
(0o |, moreover this variation together with the normal- 
ization condition for |0o) leads one to solve an eigenvalue 
problem, and a new value of (fiia)n can be calculated by 
using the resulted wave function. This should be done 
until the desired convergence of the wave-function or en- 
ergy functional is achieved. 

However, to avelliate this complication, instead of cal- 
culating the expectation value (fiia)®, we introduce a set 
of new variational parameters n^s that will play the role 
of local noninteracting occupancies (local noninteracting 
density matrices) which appear in the renormalization 
factors and double occupancies. It is then possible to let 
n ia vary independently from |0o). The energy functional 
that should be optimized has now the following form for a 
simple cubic slab geometry with periodic boundary con- 
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ditions in the x — y plane with free (001) surfaces: 
(H) = (ftaCfej + ^i)<0o|c| fc|| CT c ifc||CT |0 o ) 

(lj')fc||<T 

+ 5Z Alff (IZ^ l a ^|| ( TC J :fe|| ( T|<?!'o) - N kll nia) 

icr fcy 

+^(iV fe|| ^ n l(J - JV) + £(1 - (^ol^o)) 



The resulting non-interacting many-body ground state 
energy and wave-function are computed in the fol- 
lowing way: E NI = Y,E k , n <E P E h,n and Ivo) = 

U.B k n <E F d k h n<r4: p na\ }' where E f is the Fermi energy 
and n is the quantum number for the energy level of each 
k-point. The above non-interacting state |^o)> which is 
now implicitly a function of all the variational parame- 
ters Aicr, riicr, gi and A, should be inserted into Eq. (|10[) 
which becomes, 



where efc.. = — 2t(cosfca; + cosfcy), the Lagrange multipli- 
ers Aj<T are introduced to fix to (hi a )o- A is intro- 
duced to fix the total number of electrons, E is consid- 
ered to make sure that \ifo) is normalized, i and j are 
index of layers in the z direction and JVjfe., = N k;c Nk y is 
the total number of k-points. The optimization of the 
Lagrange function is performed through an iterative pro- 
cedure, starting with a minimization with respect to \<po), 
which leads to a Schrodinger-like eigenvalue problem that 
has to be solved for each k-point: 



J2(l^ e h + v i + ^ioOcLo-c 



liils. 



(11) 



(H) = E NI {n ia ,\ ia ,g d ,i,\ip )) - iVfc,, A i(T n lcr (12) 

i.(T 

+ A{N k]] Y n ia -N) + N h Y Ugl t n la n tlf . 



In the next step we search for the stationary points of 
the above Lagrange function of a slab geometry for a 
paramagnetic system with rii a = = m and (ni a ) = 
{^?yt)n as: 



9(H) 
dgd.i 

9(H) 
dn t 

9(H) 
9Xi 

9(H) 
dA 



2p^(U + Si^ul&Uj) + 2N k M ni 2 g hd = 0, 

Ogd,i V 9itr 

2^{U + ho±i^-h) + 2^V fe|| (A - A,) + 2N h Ug i<d 2 m = 0, 

^^Y^haCik^Wa) - 2N kll m = 0, 
fen 

(N total -2N kll Y^) = 0, 



(13) 
(14) 
(15) 

(16) 



where the spin index of renormalization factors 
and \i a is droped due to paramagnetic con- 
dition, U = £ fc || efe„<¥'o|ct fc||0 .c i jfe ||0 -|¥>o), Uj = 
_ *Sfe,| (vo IcJ fc o-Cj-fcp o- 1 ¥»o) and N to tal is the total number 
of electrons, and to.i are equal to zero at edge of 

the slab. For a detailed derivation of the saddle point 
equations for an slab geometry in the paramagnetic case 
we refer Ref. (24j. This set of nonlinear equations can 
be solved by using a nonlinear solver based on Newton 
and/or Quasi-Newton methods. Notice that \4>q) is still 
implicitly a function of the variational parameters and 



has to be updated again through Eq.( ITT]) during the 
evaluations of the saddle point equations throughout the 
optimization procedure. This means that we are all the 
time working with a |</> ) which satisfies the condition 



It should be noticed that together with the saddle point 
equations the electrostatic forces due to long-range 
electron-electron and electron-ion interactions should be 
in principle considered. However since the back ground 
permitivity of strongly correlated materials is usually 
very higb2&, we tested the solutions with various high 
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FIG. 1: (Color online) (a) Charge distribution for different 
dopings. Inset: charge transfer from bulk to surface as a 
function of doping; (b) QP weight relative to the bulk QP 
weight near the surface for different dopings. Inset shows the 
doping dependence of the bulk QP weight. 

values of background permitivity and observed that 
long-range screening is negligible^. We therefore set 
the value of the back ground permitivity to infinity in 
our calculations and ignore these effects. In order to 
numerically solve the set of saddle point equations, we 
use a 16 x 16 Monkhorst-Pack^ k-grid for which the 
total energy is well converged. We report results for 
qi as being the position dependent QP weight, which 
is a measure of the mobility of the electrons within 
Fermi liquid theory. The inverse of the QP weight is 
proportional to the mass renormalization which becomes 
divergent for qi = corresponding to an insulating 
phased. The parameters U and v are scaled by the 
tight binding parameter t. Throughout this work the 
thickness of the slab is taken L z = 90 in units of the 
lattice constant. 



III. RESULTS 

A. Hole doped correlated slab 

In Fig. [IJa) we depict the charge distribution near 
the surface for different values of doping and U = 16.2, 
which is larger than the bulk critical U for the half-filled 
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FIG. 2: (Color online) (a) Inverse correlation length as func- 
tion of doping for values of U = 16.2, 16.5, 17.0, 18.0. The 
inset shows the inverse of the correlation length as function 
of U for four values of 5=0.009, 0.013, 0.017, 0.032 from but- 
torn to top curves. 



jbuik,hf _ -^g rpj ie sur f ace re gi on i n which 

the charge density recovers its bulk value is doping de- 
pendent, resulting in the doping dependent correlation 
length. Higher doping corresponds to lower correlation 
length. In the inset of FigOJa) we present the charge 
transfer from the bulk to the surface (n sur f ace — ribuik)- 
The doping dependence of this charge transfer is non- 
monotonic and is maximum around S = 0.15. While 
our results for the charge transfer are in agreement with 
recent DMFT calculations for a hole doped semi-infinite 
single band Hubbard models in the limit of large enough 
doping, our scaling analysis shows that considering only 
few layers for the QP calculation may not be enough, spe- 
cially for values of doping near half filling for which the 
correlation length is larger that 6 lattice constants. In 
Fig. [ljb) the spatial distribution QP weights (qi — qbuik) 
are plotted for different values of doping and U = 16.2. 
Like in the half-filled caseZ&21 the QP of electrons near 
the surface sites is suppressed due to the reduced coordi- 
nation number together with the charge transfer to the 
surface sites from the bulk, which in turn results in a 
lack of kinetic energy and an enhancement of correlation 
effects. One can also observe Friedel oscillations which 
are more pronounced for higher doping due to lower cor- 
relation lengths. The inset of Fig. [IJb) shows the dop- 
ing dependence of the bulk QP weight, qbuik, which is 
in agreement with previous works and shows that by in- 
creasing the doping, correlation effects are weaker—. The 
correlation length can also be extracted from the spatial 
distribution of the QP weight near the surface. Similar 
to the dependence of the charge density, the QP weight 
recovers its bulk value within a characteristic length scale 
that depends on the correlation length. Friedel oscilla- 
tions can also be observed but are suppressed for lower 
doping. Following Ref. we observe that the spatial dis- 
tribution of \J q(x) — y/qbuik is we ll fitted by an exponen- 
tial decay for different values of the Hubbard repulsion 
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FIG. 3: (Color online) QP weight distribution and charge 
distribution (inset) for U = 16.2, S — 0.002 and three different 
values of the electric field. 



c 
Q 

3 

TD 
>» 
CO 

c 

CD 

Q 

a> 
E? 

O 




20 30 40 50 60 70 
Index of atoms in z direction 



90 



FIG. 4: (Color online) Charge distribution; the inset shows 
the charge transfer as function of U for fixed v = 0.88 and 
5 = 0.002 



and doping: 

Vl(x) = y/q bu lk + {y/q S urface ~ \/%ulk)e~^ X ~ l \ (17) 

where £ is the correlation length and x the number of 
layer, starting from x = 1. Since the correlation length, 
£, depends on both U and S, by fitting separately the 
spatial distribution of the QP weight we extract the cor- 
responding correlation lengths. The results are summa- 
rized in Fig. [2j where l/£ is plotted as a function of dop- 
ing for different values of the Hubbard repulsion. We 
can extract a simple power-law dependence for the in- 
verse correlation length: | = AS V , with a mean-field-like 
exponent&2£, 77 = 0.5 ± 0.07, and a prefactor A that is 
only a function of U . The inset of Fig.[2]shows the inverse 
correlation length as a function of U for different dopings 
and, as expected, it is enhanced for higher Hubbard re- 
pulsions. The power law dependence of the correlation 
length versus doping shows that for half-filling the cor- 
relation length diverges which is a signature of the MIT 
that occurs for U > U^ ulk - h f . A similar dependence of 
the correlation length versus Hubbard repulsion is ob- 
served for half-filling but when U < Jjbuik,hf25_ j n t j ie 
latter case the criticality is governed by the Hubbard re- 
pulsion rather than the doping level. 

B. The effect of electric field 

The effect of an external electric field perpendicular to 
the slab on the spatial distribution of the QP weight is 
shown in Fig.|3]for U = 16.2, 5 = 0.002 and different val- 
ues of the voltage difference. The inset shows the charge 
distribution for the same parameters. The main effect 
of the electric field is to redistribute the charges within 
the slab, however in the strongly correlated regime when 
the Hubbard repulsion exceeds a certain crossover value, 
correlation effects enhance the charge transfer from less 
correlated sites to more correlated ones. This correlation 



enhanced charge redistribution results in the accumula- 
tion of charges near the surface layers, bringing one side 
of the slab very close to half-filling. This effect is largest 
for U > Uj! ulk ' hf . 

To better clarify the correlation effects on the surface 
states of a correlated slab in the presence of an electric 
field we depict in Fig.@]the charge and quasiparticle(QP) 
distribution of a slab consisting of 90 layers thick and a 
voltage difference v = 0.88. The charge distribution for 
U = 15.22 shows peaks near the surfaces, as expected, 
however this behavior disappears for U = 15.74 and 
U = 16.2. This shows a clear crossover regime related 
to the enhancement of correlation effects. On the other 
hand the naive expectation that the effect of an increased 
Hubbard repulsion is only to screen out the electric field, 
fails to explain the behaviour of the system in the pres- 
ence of the electric field in the strong coupling regime. 
As shown in Fig. [4] by increasing the Hubbard repul- 
sion, the charges do not go away from the surface but 
instead are accumulated at the surface. This mechanism 
of charge transfer from the places with enhanced der- 
ealization to the places with enhanced correlations leads 
to a non-trivial enhancement of QP difference between 
the surfaces for large Hubbard repulsions. To further un- 
derstand the charge redistribution enhancement due to 
correlation effects, we present in the inset of Fig. 2] the 
charge difference between the layers with highest charge 
density and the layers with lowest charge density as func- 
tion of U. This can be considered as a measure of the 
charge transfer throughout the system. As is clear from 
the inset of Fig. @] there is a crossover value for U, given 
a fixed doping S = 0.002 and voltage difference v = 0.88. 
Above this value the effect of the U plays a different role 
in the charge redistribution in the system. While below 
the crossover interaction the Hubbard repulsion competes 
with v to prevent charge redistribution due to voltage 
difference, above the crossover it enhances the charge re- 
distribution in favor of v. As is obvious from Fig. [3] the 
maximum QP weight is already achieved after a few lay- 
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FIG. 5: (Color online) QP for x=l and x=90 as function of 
voltage difference for U — 16.2 and three values of 5. 

ers from the surface on that side of the slab where the de- 
viation of the charge density from half-filling is maximal. 
The reason that the QP weight is not maximal exactly 
at the surface is because of the suppression of the kinetic 
energy near the surface. On the other side of the slab, for 
larger electric fields the charge transfer assures that the 
charge density is near half-filling. Therefore, due to local 
correlation effects the QP weight is strongly suppressed. 
While the charge density near the surface is very close to 
half-filling (i. e. n — 1 ~ 10~ 7 ) one may infer that the 
residual QP indicated in Fig. [3]for x = 1 is mostly due to 
the proximity of the surface site to sites with higher QP 
weight rather than due to the local doping effect of these 
regions. In order to better understand the dependence 
of the QP weight on opposite sides of the slab, in Fig. [5] 
we show the QP weight for layers x = 1 and x = 90 as 
function of voltage difference for three different values of 
doping. 

The QP weights on the two surfaces differ by many 
orders of magnitude. For larger doping, higher electric 
fields are needed in order to achieve the same QP weight 



difference. This is because of the competing influence of 
doping and Hubbard repulsion on the correlation effects. 
The huge difference in QP weight near the two surfaces 
could be used for creating a transistor-like device made of 
strongly correlated materials. By using the surface states 
to conduct current one can simply switch on/off the de- 
vice by switching the polarity of the gate. Thus, turning 
on/off the electric conduction is now a consequence of 
the surface resistance ratio of the two sides. 



IV. CONCLUSIONS 

By using an inhomogeneous Gutzwiller approach ap- 
plied to the paramagnetic single band Hubbard model 
for a slab geometry we described a hole doped Mott thin- 
film. In the absence of applied electric field we calculated 
the position dependent charge density and QP weight and 
showed that the inverse correlation length has a power 
law dependence on doping. 

When a perpendicular electric field is applied, charges 
will accumulate on one side of the slab. This correla- 
tion enhanced charge redistribution will in turn induce 
a large difference in the QP weight on the two sides of 
the slab, which was found to be as large as five orders 
of magnitude. We propose that a three terminal device 
with surface contacts can take advantage of this effect. 
For resistance switching purposes one would expect large 
on/off ratios of surface resistances when the electric field 
switches polarity. 
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